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Abstract 

We derive the continuum equations and boundary conditions governing phonon- 
mediated heat transfer in the limit of small but finite mean free path from asymp¬ 
totic solution of the linearized Boltzmann equation in the relaxation time approx¬ 
imation. Our approach uses the ratio of the mean free path to the characteristic 
system lengtliscale, also known as the Knudsen number, as the expansion parameter 
to study the effects of boundaries on the breakdown of the Fourier descrition. We 
show that, in the bulk, the traditional heat conduction equation using Fourier’s law 
as a constitutive relation is valid at least up to second order in the Knudsen number 
for steady problems and first order for time-dependent problems. However, this de¬ 
scription does not hold within distances on the order of a few mean free paths from 
the boundary; this breakdown is a result of kinetic effects that are always present in 
the boundary vicinity and require solution of a Boltzmann boundary-layer problem 
to be determined. Matching the inner, boundary layer, solution to the outer, bulk, 
solution yields boundary conditions for the Fourier description as well as additive 
corrections in the form of universal kinetic boundary layers; both are found to be 
proportional to the bulk-solution gradients at the boundary and parametrized by the 
material model and the phonon-boundary interaction model (Boltzmann boundary 
condition). Our derivation shows that the traditional no-jump boundary condition 
for prescribed temperature boundaries and no-flux boundary condition for diffusely 
reflecting boundaries are appropriate only to zeroth order in the Knudsen number; 
at higher order, boundary conditions are of the jump type. We illustrate the util¬ 
ity of the asymptotic solution procedure by demonstrating that it can be used to 
predict the Kapitza resistance (and temperature jump) associated with an interface 
between two materials. All results are validated via comparisons with low-variance 
deviational Monte Carlo simulations. 
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1 Introduction 


Microscale and nanoscale solid state heat transfer as mediated by phonon transport has 
received considerable attention in connection with a number of diverse practical applica¬ 
tions, such as heat management in microelectronic devices, passive cooling and thermo¬ 
electric energy conversion 11, but also due to the number of scientific challenges it poses. 
Particularly notable is the wide range of scales present in these problems, typically start¬ 
ing from the atomistic (including quantum) and extending to the macroscopic (device). 
Kinetic-theory approaches based on the Boltzmann transport equation (BTE) [2|, espe¬ 
cially if informed by ab-initio information on the material properties [3jj5], can be quite 
effective in bridging this range of scales. One limitation of such approaches appears in 
the small mean free path limit, (Kn) <C 1, where kinetic descriptions become stiff. Here, 
(Kn) denotes the Knudsen number defined as the ratio of the mean free path to the 
characteristic system lengthscale; a more precise definition will be given in section [2] 

As is well known, in the limit (Kn) —> 0, the stiff Boltzmann description need not be 
used because it can be replaced by the heat conduction equation; derivation of the bulk 
thermal conductivity from the Boltzmann equation in the relaxation approximation 
via a Chapman-Enskog type of expansion 16 ,[T] is well established, thus providing a 
“pathway” for recording the effect of molecular structure on the constitutive behavior 
in that limit. However, the Chapman-Enskog expansion is only applicable in the bulk 
and provides no information on the boundary conditions that need to supplement the 
heat conduction description in order to obtain solutions that are consistent with the 
(more fundamental) Boltzmann solution. Moreover, a rather large gap exists between 
lengthscales that truly satisfy (Kn) —>• 0 and the regime where Boltzmann equation 
solution is no longer problematic ((Kn) > 0.1). 

In this paper, we use an asymptotic expansion procedure using (Kn) as a small 
parameter to derive, from the BTE, the “continuum” equations governing phonon- 
mediated heat transfer in the small mean free path limit. This procedure recovers the 
classic heat conduction equation (including Fourier’s law as a constitutive relation) as 
the equation governing the temperature field that is consistent with solution of the 
Boltzmann equation to order (Kn)°, as expected. However, in contrast to Chapman- 
Enskog-type procedures, this procedure, also derives the boundary conditions that the 
heat equation is to be solved subject to. Specifically, for fixed temperature boundaries, 
the Fourier boundary conditions are found to be of the Dirichlet type at the boundary 
temperature; for diffusely specular walls, the Fourier boundary conditions are shown 
to be the Neumann no-flux boundary condition. Although these results have been 
empirically established centuries ago, this is the first time they are shown to arise, 
rigorously, from a solution of the Boltzmann equation. 

More importantly, by extending the asymptotic expansion to first and second order 
in (Kn), we derive the governing “continuum-level” equation and boundary conditions 
for finite but small values of the Knudsen number ((Kn) <C 1). Specifically, for steady 
problems, the governing equation is shown to be the steady heat conduction equation 
up to order (Kn) 2 , while the corresponding boundary conditions are shown to be of the 
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temperature-jump type, with jump coefficients that, in general, depend on the material 
and boundary properties. For unsteady problems, we show that the governing equation 
is the unsteady heat conduction equation up to first order in (Kn) with boundary 
conditions remaining the same as in the steady case up to that order for the case of 
prescribed-temperature boundaries. 

Jump boundary conditions have been observed before in solutions of the Boltz¬ 
mann equation |8y9| and attempts were made [8] to explain these invoking differences 
in local equilibrium conditions across interfaces. The present work shows how tem¬ 
perature jumps arise as a result of the incompatibility between the isotropic distribu¬ 
tions associated with boundary conditions and the anisotropic distribution associated 
with non-equilibrium resulting from transport (temperature gradients). A well-known 
manifestation of this physical behavior are the temperature jumps associated with the 
Kapitza interface problem. In section [8] we show how our asymptotic approach can be 
used to calculate the interface conductance (and associated temperature jump) from first 
principles (at the kinetic level, that is, given the interface transmission and reflection 
coefficient). 

The temperature jump relations derived in this work are manifestations of what 
is known in the kinetic theory community as “slip”, which gives its name to the slip 
regime, 0 < (Kn) <0.1. It is generally known 10 UK that in this regime the material 
constitutive law may still be used unmodified and kinetic effects are accounted for by 
modified boundary conditions. In the field of rarefied gas dynamics, Cercignani |12| 
and Sone with co-workers 13,[l4j were the first to provide systematic asymptotic solu¬ 
tions up to second order in (Kn), demonstrating the possibility of using the traditional 
“continuum” fluid dynamics, albeit with modified boundary conditions, beyond the slip 
regime and into the early transition regime. The transition regime is typically defined 
by 0.1 < (Kn) < 10 and represents the regime in which transport transitions from diffu¬ 
sive ((Kn) <C 1) to ballistic ((Kn) S> 1). Discussions of the use of asymptotic solutions 
of the Boltzmann equation in rarefied gas dynamics can be found in 1 10,15,16 1 . 

The practical implications of the present work are twofold: first, solution of the heat 
equation is significantly easier (analytically or numerically) compared to the Boltzmann 
equation, especially in the regime (Kn) -C 1 where the latter becomes stiff. In addition 
to ease of solution, centuries of investment in continumm formulations such as the 
heat equation, either in the form of education, mathematical solution techniques or 
numerical solution software, make this by far the preferred approach. This can be 
easily seen from the considerable efforts expended in developing approximate ’’effective 
thermal conductivity” concepts that enable the use of Fourier’s law in the transition 
regime. The present work provides rigorous methods for obtaining solutions consistent 
with the Boltzmann equation in the slip and early transition regime. Studies in rarefied 
gas dynamics show that, depending on the problem and the amount of error that can 
be tolerated, slip/jump formulations could be used up to (Kn) ~ 0.5 and sometimes 
beyond [lT]. Second, by using the asymptotic solution as a control in deviational Monte 
Carlo schemes, one can overcome the stiffness associated with the (Kn) <C 1 regime. This 
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happens because 18,19 the asymptotic solution becomes increasingly more accurate 
as (Kn) —> 0, thus requiring increasingly less computational resources to describe the 
deviation therefrom as this limit is approached. This yields computational methods 
that are able to efficiently simulate problems characterized by (Kn) < 0.1 locally or 
globally, in contrast to traditional Boltzmann solution methods. 

The present paper is organized as follows: in section [2] we introduce the governing 
(Boltzmann) equation and the notation used in this paper; in section [3j we present the 
asymptotic analysis leading to derivation of the governing equation in the bulk up to 
second order in the Knudsen number. Associated boundary conditions and boundary 
layer corrections up to first order in the Knudsen number are derived in section [4j In 
section [5] we present results obtained from extending the boundary layer analysis to 
second order in Knudsen number. In section [6] we summarize and discuss our results 
and provide example applications to one-dimensional and two-dimensional problems. In 
section [7] we discuss the applicability of the asymptotic theory and its results (govern¬ 
ing equations, boundary conditions and corrective boundary layers) to time-dependent 
problems. In section [8] we show how the asymptotic theory can be used to calculate 
the Kapitza conductance (and temperature jump) associated with the interface between 
two materials. We conclude with some final remarks in section [9j 


2 Background 

We consider the Boltzmann equation for phonon transport in the relaxation time ap¬ 
proximation 

d f f loc - f 

t ^+ V ,- V x ,/=4 -( 1 ) 

at' t(lu,p,T) 

where / = f(x.',uj,p,tt,t') is the occupation number of the phonon states, x' the po¬ 
sition vector in physical space, V g (u,p) the group velocity, u the phonon frequency, 
p the phonon polarization, Q the unit vector denoting phonon traveling direction, T 
the temperature and / loc an equilibrium distribution at the “pseudotemperature” T\ oc 
defined by energy conservation considerations (refer for instance to [8,201 for details on 
the definition of / loc ). 

In this work we primarily consider steady problems. Extension to time-dependent 
problems directly follows by extending the methodology presented here. Scaling analysis 
in section [7] shows that, assuming diffusive time scaling, time dependence may modify 
the results presented here at order (Kn) 2 . In other words, the results obtained for 
steady state in this paper may be applied directly to order (Kn)° and (Kn) 1 with very 
few modifications, explained in section [7j 

Assuming small deviations from equilibrium at temperature T eq , the linearized steady- 
state Boltzmann equation reads 


v 9 -v x ,/ 


d 


£(/ d ) - f d 

r(a;,p,T eq ) 


( 2 ) 
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where f d = (/ - / eq ), with / eq = [exp(hu / k h T eq ) - 1] 1 . 

By noting that £(/ d ) = (T[ oc — T eq )df eq /dT and writing energy conservation 
the form 

f cu*)™*/-! 

J W f ,p' j id',p' 


20 


m 


r T 4vr 

where D = D(uj,p ) denotes the density of states, we obtain the expression 

d ,_Ly,n'^£dWdu;'dr« 


(3) 


£(/ d ) = 


a 


dT 


(4) 


Here, and in what follows, unless otherwise stated, r 
expression, 


C T 


f Dhw df eq , 

L P T dT W 


t(uj,p, T eq ). In the above 


(5) 


Also, 11 and d 2 H respectively refer to the unit vector defining the direction of propaga¬ 
tion and to the differential solid angle, expressed as sin (d)dOd(j) in spherical coordinates. 
In the interest of simplicity, in the above expressions and in what follows, we use a 
single integral symbol to denote both integrals over multiple variables and sum over 
polarization. 

In this study, relaxation times and group velocities may depend on frequency and 
polarization. For this reason, the Knudsen number is defined in an average sense. We 
choose the following (somewhat arbitrary) definition 


(Kn> 



( 6 ) 


where 

C U , P = (7) 

and Kn w .p = A U ,p/L = V g (uj,p)T(uj,p,T eq )/L, which we will denote by Kn. In the 
expression for Kn^p, V g (uj,p ) = ||V 9 (w,p)|| is the magnitude of the group velocity 


3 Asymptotic analysis for the bulk 

Introducing the dimensionless coordinate x = x'/L as well as the normalization 

f 


<F = 


cp q 

dT 


( 8 ) 


we write the Boltzmann equation in the form 

£($) - $ 


n ■ v x <h = 


Kn 


(9) 
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where 


f „ ^ <S>d 2 Clduj 

47TT 


£($) = 


c T 


( 10 ) 


The usual macroscopic quantities of interest such as temperature, energy density 
and heat flux can be calculated from 


Ttot — T ea + 


1 


1 tot — ^eq t . ^ I v - y ^4> 

47r ° Ju;,p,n 


Cu,p&d 2 £lcbj = T eq + T(x) 


LU,P,tt 


C UtP &d 2 flduj 


E t °t - E eq + 

q" = — [ C u p V„^Ctd 2 nduj 
4vr 7 WiPi n 


( 11 ) 

( 12 ) 

(13) 


We will refer to T(x) as the deviational temperature, since it represents the deviation 
from the equilibrium temperature T eq . 


3.1 Bulk solution 


The asymptotic solution relies on a “Hilbert-type” [21] expansion of the solution d> in 
the form 

OO 

T = ^(Kn)"$ n (14) 

n =0 

Given the nature of the proposed solution, similar expansions can be written for the 
temperature and the heat flux fields 


T = En=o <Kn)"T n 

q" =£n=o<K n > n q" 


(15) 


In this section, we only consider solutions far from any boundary. As will be shown 
below, close to the boundary, kinetic effects become important due to the incompatibility 
of the bulk solution with the kinetic (Boltzmann) boundary condition and a separate, 
boundary layer analysis is required. Therefore, we let = ^(Kn) n 4>G n be the bulk 
solution, anticipating that where represents kinetic boundary layer 

corrections that are zero in the bulk and will be similarly expanded later. When the 
expansion for 4 >g' is inserted in the Boltzmann equation we obtain 


ft • V x f](Kn)"$ Gn = f](Kn) n ^ (16) 

n= 0 n= 0 

By equating terms of the same order ((Kn) 1 and higher powers) and assuming that 
Kn ~ (Kn), we obtain the following relationship for all n > 0 

B • V x <l>Gn = [^(^Gn+l) — ^Gn+l] ■ (17) 
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In addition, considering the two terms of order 0 in the right hand side of (16), we find 
that <h G o is determined by the solution of the equation 


L P ,n G0 d 2 flduj 


^GO = £(<hco) = 




cv 


(18) 


The assumption Kn ~ (Kn) is easily satisfied when the range of free paths is rel¬ 
atively small (and is exactly satisfied in the single free path case A^ p = A = const), 
but becomes harder to justify in materials with wide range of free paths. In the latter 
cases, it has the effect of reducing the value of (Kn) for which the theory presented here 


is valid. This is further discussed and quantified in section 4.1.1 
From equation 


we deduce that <l> G o is a function that depends on x only, since 
this is the case for £(4 >g , o)- We note here that any function that only depends on x is 
a solution. Additionally, since <F G o = 4> G o(x), we find that the zeroth order deviational 
bulk temperature is given by 


Tco(x) = —^ [ C^p^Go(^)d 2 ftduj = $go(x). (19) 

and that 

Qgo = 7- [ C^pK g T G o(x)nd 2 fidu; = 0 (20) 

47r Juj,p,n 

At this stage, the spatial dependence of <h G o is undetermined. The additional infor¬ 
mation needed will be inferred from the application of a solvability condition to <h G i. 
Using © we find the following expression for the order 1 solution 

Kn 

$G1 = £($Gl) ~ ixr ■ ^X^GO (21) 

(Kn) 

This equation states that a necessary condition for <F G i to be the order 1 solution is 
that it is equal to the sum of —Kn(Kn) _1 fi • Vx^go an d a function that only depends 
on x. Since the temperature associated with ft ■ V x 4>go is zero, we can write 


<f>Gi = Tg\ ~ Kn(Kn) 1 ft ■ V x Tq 0 (22) 

Finally, order 2 may be derived following the same procedure for eq © for n = 1, 
which yields 

Kn Kn 2 

$ G2 = £(4> G2 ) - — ft ■ V x r G1 + —2 n • v x (ft ■ V x T G0 ) (23) 

(Kn) (Kn) z 

In the following section, while deriving the governing equation for Tqq, we also show 
that the temperature associated with <I> G2 is £(4> G2 ) = T G2 . 
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3.2 Governing equation for the temperature field 

The solvability condition required to determine 4>Gn is the statement of energy conser¬ 
vation ([3]) which, applied to ^Gn+n becomes 


C, 


Ld,P 


f UJ,p 


£{<5>)du = [ dud 2 fl 

Jw, P ,n 4ttt 


(24) 


Using ( |17[ ) results in the following condition 

/ C U:P V g fl ■ V x ^ Gn dud 2 fl = 0 (25) 

J co,p,n 

that needs to be satisfied for all n > 0. Applying this relationship to <3?gi, we obtain 




Kn 

Cuj^Vgfl ■ V x ( Tqi — fS • V x 7go ) du;d z fl = 0 


(Kn 


(26) 


In the above expression, the integral over the solid angle is zero in all terms where a 
component of the traveling direction appears with an odd exponent. This implies 

V 2 T G0 = 0 (27) 

This concludes the proof that the 0-th order temperature field obeys the steady state 
heat equation. Moreover, from (23) it follows that 

Kn 2 


$G2 = Tg-2 - 


Kn 

W) 


• V x 7gi + 


(Kn) 


fi • V x (f2 • V x Tgo) 


(28) 


In Appendix [A] we show that higher-order (in fact, possibly all order) terms similarly 
obey the heat equation. In other words, Xgi(x) and 7g2(x), are determined by solution 
of 


VxTgi = 0, V^Tg 2 = 0 


(29) 


Before we close this section, we note that although in the Laplace-type equations 
derived above for the temperature the thermal conductivity does not appear, the above 
asymptotic analysis still clearly predicts that in the bulk, the material constitutive re¬ 
lation (thermal conductivity) is equal to the ’’traditional’' bulk value. This can be seen 
from first-principles by inserting (22) into (13) to obtain 


(Kn)qGi = -7- 


V 2 t 

v 9 1 


.1 oj, P) Cl L 


Cuj t p£l (Li • V X T 0 ) duicP^l = — kV x 'Tgo 


(30) 


where the second equality follows from recognizing the well known expression 

1 


K = 


VgTC U:P du; 


> LU,p 


(31) 











4 Order 1 boundary layer analysis 


In this section, we extend the asymptotic analysis of the previous section to the vicinity 
of boundaries, where as will be shown below, a boundary layer analysis is required 
for matching the bulk solution of the previous section to the kinetic (BTE) boundary 
conditions of interest. Here we will consider two kinetic boundary conditions, namely, 
those of prescribed temperature and diffuse adiabatic reflection. In this work we assume 
that boundaries are flat; boundary curvature will be considered in a future publication. 
Without loss of generality we assume that the boundary is located at x\ = 0 and with 
an inward normal pointing in the positive x\ direction; X 2 and X 3 will denote cartesian 
coordinates in the plane of the boundary. Moreover, we will use Hi, and H 3 to refer 
to the components of the unit vector H in the coordinate system (xi,X 2 ,xs). In other 
words, Hi = cos( 0 ), H 2 = sin( 0 ) cos( 0 ) and H 3 = sin( 0 ) sin( 0 ). 

We now derive the general equation governing the boundary layer correction required 
in the boundary vicinity for matching the bulk solution to the kinetic (BTE) bound¬ 
ary conditions. We introduce the boundary layer function <&k- written as a Hilbert 
expansion (<h n = &Gn + &K 71 ) with d>A'o = 0 and insert it in the Boltzmann equation, 
obtaining 


00 3 


yy<Kn)% 


d$Ki 

dxj 


E< Kn>‘ 
2=1 


C{^Ki) ~ ®Ki 
Kn 


(32) 


In the vicinity of the boundary, a new characteristic lengthscale, namely the distance 
from the boundary, becomes important. Similarly to [lb] , we introduce a “stretched” 
variable defined by rj = xi/(Kn). Equation (32) can thus be written in the form 


2=1 


d§i<i 

dr] 


E< Ki 

2—1 


y C-i^Ki) - $Ki 
' Kn 


E< Ki 

2=1 


OX 2 OX 3 J 


( 33 ) 

By equating terms of the same order, we find that each boundary layer term is solution 
to a ID (in physical space) Boltzmann-type equation. For ‘h^i) this equation is 


Hi 


thl 1 k 1 
dr] 


(Kn) 


Kn 


(34) 


The equations for &Kn, n > 2 include ’’volumetric source” terms resulting from the 
derivatives of the lower order boundary layers in the boundary tangential directions (X 2 
and X3). Specihcally, for each order i > 2: 


n d$Ki ~ &Ki d<f>Ki -1 n 

n '^T = <Kn> -Kn- +Si 3 ^S-j 


(35) 


The case i = 2 will be considered in the following section, where second-order boundary 
layer analysis is carried out. 
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4.1 Boundary conditions for prescribed temperature boundaries 


The term “prescribed temperature boundary” is typically used to describe a boundary 
approximating a black-body, absorbing incoming phonons and emitting phonons from 
an equilibrium (isotropic) distribution at a given temperature. In other words, the 
Boltzmann boundary condition associated with such a boundary at deviational temper¬ 
ature T b is a Bose-Einstein (equilibrium) distribution at the wall temperature, denoted 
here by / eq (o;; T eq + T b ). In the linearized case, the incoming distribution of deviational 
particles is therefore 

df eq 


fb = T b 


dT 


(36) 


or simply, in terms of quantity 4> defined in ([ 8 ]), d>t, = T 5 . 

We note that 4 >go is isotropic and is thus able to match 4>{, provided we set Tq 0 = T b 
at the boundary. Therefore, at order 0, the solution to the Boltzmann equation with 
prescribed temperature boundaries is given by the heat equation complemented by the 
traditional Dirichlet boundary conditions and no boundary layer correction is required 
(dWo = 0 , which also implies that To = Tq 0 ). 

This situation changes at order 1. The order 1 distribution 4 >gi = Tq\— Kn(Kn) _ 1 fi- 
V x Tgo is not isotropic due to the gradient of Tq o- As a consequence, there is a mismatch 
between the order 1 solution and the boundary condition (which has been satisfied 
by 4 *go and is thus zero for all subsequent orders). This mismatch can be corrected 


by introducing a boundary layer term 4 >a^i governed by equation (34) and subject to 


boundary condition 4 >_r'i|??=o + 4>gi|?j=o = 0, which translates into the following relation 


4 >a'i 177=0 = —Tgi\ v =o + 


Kn 

(Kn) 


B ■ VxTgoI^o 


(37) 


The term V x Tgo|? 7 =o is known from the order 0 solution. The term Tgi\ v =o is 
unknown and determined by the fact that there exists only one value for Tq\ L = q such 
that 4>a;i tends to 0 for 7/ — > 00 1 101. This determination proceeds by writing 4>au = 
4>a'i.i + 1,2 + &ki ,3 where each of &Ki,i, z = 1 , 2 ,3 is the solution to an equation of 


the form (34) with the associated boundary condition: 


4 * A'l I n =0 = ~°i + 


— fl 
(Kn) 1 


dT G0 


dx; 


(38) 


r ]=0 


Anticipating the values of 4>A'i,i to scale with dTGo/dxi\ v= o in the above equations 
we have set Tgi\ v =o = Yli c i(dTGo/dxi)\ v =o- The constants ci, C 2 , C 3 are uniquely 
determined by the condition that 4>a'i,i, 4>ai,2 and 4>a _ i,3 individually tend to zero for 
rj —> 00. 

One can easily verify that for i = 2,3, Cj = 0, with 


4 > AU,j = 4- 


Kl,i 


rjrp I Kn n. 9Tgo 

01 go _ I (Kn)“ l dxi 


dXi 


rj=0 


exp 


(Kn)77 

KnQi 


0 , 


for > 0 
for Oi < 0 


(39) 
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is a solution to (34) with boundary condition (38). The temperature field associated 
with these functions is zero. Here we note that the above solutions have the property 
£(^ 1 , 2 ) = £(^' 1 , 3 ) = 0 and thus are also solutions of (34) with the term C(&ki) 
removed. We will use this observation throughout this paper for obtaining analytical 
solutions to a number of boundary layer problems. 

The problem for 4>xi,i must be solved numerically. Given the boundary condi¬ 
tion it needs to satisfy, we write ( ^Tgo/ dx\)\ Q and solve for ^ki,i- 

The numerical method developed and used for this purpose is explained in detail in 
Ref 22 . In the case of a Debye and gray material referred to here as the single 
free path case (Kn = (Kn) for all u,p), it yields c\ = 0.7104, while the resulting 
tk 1,1 = J u n Ki^ducPSl/AiT is plotted in Figured! We note that Refs 
also report the value 0.7104 in the context of other kineuc particle transport 
velop other efficient methods for solving this problem. 


23 


and 


24 

He 



Amax 


Figure 1: Temperature profile associated with tki,i = £ru,i/(<£7go/<9xi)|?7=o, for three 
relaxation time models. The x-axis is scaled by the maximum free path A max of each 
model. 


In summary, the boundary condition for the order 1 bulk temperature held is 


rp ( n s 9T G o 

Tc ' {x ‘ = 0) = C1 fcT 


(40) 


or more generally 


Xb 

where dTco/dn refers to the derivative in the direction of the normal to the boundary 
pointing into the material, n, and x& the boundary location. In other words, the bound¬ 
ary condition is of the jump type and the associated temperature jump is proportional 
to the derivative of the 0th order solution in the direction normal to the boundary. 



£(311 x= ci 


dT G q 
dn 
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The amplitude of the corrective boundary layer that is added near the wall is also 
proportional to the normal derivative: 


Tki,i = tki, i 


dT G0 


dn 


(42) 




Note that although a non-zero temperature field is associated with &ki, i ; the corre¬ 
sponding heat flux is zero. This is explained by the fact that < 3 ?_r'i,i ) by construction, 
tends to 0 at infinity. Since the boundary layer problem is one-dimensional in space, 
by energy conservation, the heat flux has to be constant in x\ and is therefore zero 
everywhere. We also note that although and 4^1,3 do not contribute to the tem¬ 

perature field, they do contribute in the heat flux in the direction parallel to the 


boundary. Their contribution can be obtained by substituting (39) into (13); the result 
is summarized in section 1 


4.1.1 Numerical solution for complex material models 


In section 4.1 we reported the value of the coefficient c\ and boundary-layer function 


in the single free path case. In this section we report results for two more realistic 
material models. Specifically, we consider a material with realistic dispersion relation 
and a single relaxation time, as well as a material with realistic dispersion relation and 
frequency-dependent relaxation times. The dispersion relation in both cases is taken 
to be that of the [100] direction in silicon. The single relaxation time is taken to be 
40ps. In the case of a variable relaxation time we use a slightly modified Born-von 
Karman-Slack (mBvKS) model [25j| with parameters from 126 and 118], where the grain 
size used for boundary scattering is 0.27 mm instead of 2.7 mm. The reason for this 
approximation is that it facilitates the verification of the order 1 behavior with Monte 
Carlo simulation. We do not consider optical phonons in this work, but the method can 
be straightforwardly extended to this case. 

We find c\ = 1.13 in the single relaxation time model and c± = 32.4 in the mBvKS 
model. The associated boundary layers are plotted in figure [l] It is important to note 
that: 


- The values of coefficient ci and the function tk i,i depends on the definition of (Kn) 
or, equivalently, (A), which is rather arbitrary. This, however, does not influence 
the final result because the asymptotic temperature field, ultimately (see (|15[)) 
depends on the products ci(Kn) and rxi,i(Kn) (see for instance solution (|88|) ). 


- The boundary layer in the mBvKS model is particularly wide (on the order of 
millimeters). This observation, as well as the large value of ci, is a manifestation 
of the stiffness (multiscale nature) of this problem, resulting from the wide range of 
free paths present in this material; mathematically, it is due to the factor Kn/ (Kn) 
that appears in (37) and which tends to give more weight to modes with very 
large free paths and makes the assumption Kn ~ (Kn) hard to satisfy. Since, by 
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assumption, the sum of all fg„(Kn} n should exist -which requires <L n (Kn) n < 1- 
this has the overall effect of limiting the range of applicability of the asymptotic 
model to Knudsen numbers that are lower than the nominal (Kn) < 0.1. It is 
important to note, however, that this limitation is a result of the fundamental 
physics of the problem: even at “low” Knudsen numbers given by (Kn) < 1/ci, 
there exist modes with long free paths (i.e. Kn ~ 0(0.1)) introducing kinetic 
effects and making the zeroth order solution (V^Tgo = 0 ) inadequate. 

4.1.2 Validation 

We validate our result using a one-dimensional problem, in which a mBvKS material 
is placed between two boundaries at prescribed temperatures and located at x\ = — L 
and x\ = L, respectively. The order 0 (traditional Fourier) solution to this problem is a 
linear temperature profile Tq o(*i) which yields a heat flux ksi-m^Tgo/L, where AX^o 
is the temperature difference between the boundaries; here, ksi_m denotes the bulk 
thermal conductivity associated with the mBvKS material. The temperature profile 
Tgi is obtained by solving the Laplace equation with jump conditions 

T g i(xi = =Fl) = ±ci ^ (43) 

and yields the modified heat flux Ksi-M(l — ci(Kn))AXGo/X. We note that when cal¬ 
culated from an order n temperature field, the heat flux is inherently an order n + 1 
quantity; in other words, the above result is correct to order 2. In Figure [2j we plot the 
difference between the actual heat flux (q". , obtained using deviational Monte Carlo sim¬ 
ulation jl8j27]) and the asymptotic approximation, both normalized by ksi-mATgo/L, 
namely, e q = q'^L /(ksumATgo) — (1 — ci(Kn)). The observed asymptotic behavior is 
order 2 which validates the order 1 accuracy of the asymptotic solution. 


4.2 Boundary condition for a diffuse adiabatic boundary 

The case of diffuse adiabatic boundaries can be treated through a similar approach, 
where the mismatch between the bulk asymptotic solution and the boundary condition 
is analyzed and corrected. The boundary condition at the kinetic level is given by (281 

<h|x 6 = -- [ n\d 2 n’ for Lb > 0 (44) 


A major difference from the prescribed temperature boundary is that applying this 
condition to the Oth order bulk solution gives no information, because <3 ?go satisfies 
(44) regardless of its value at the wall. The boundary condition for Tqo is obtained by 
analyzing the order 1 mismatch. The order 1 boundary layer problem may be defined 
by applying the boundary condition (44) to $1 = 4>gi + $ki- It results in the following 
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Figure 2: Validation of the first order asymptotic theory for prescribed temperature 
boundaries. The solid line denotes the normalized (by the Oth order, traditional Fourier, 
result) difference between the heat flux predicted by the asymptotic theory and MC 
simulation results. The dashed line denotes a slope of 2. 


condition: 


Tg\ |t7=0 - ft ■ V x TGo|r;=0 + 4>A'l|?7=0 = 

(Tgi\ v =o — ^ ■ VxTgoI^q + ^/n|?7=o) for > 0 


(45) 


77 Jn ’ ± <o 

The isotropic term T G \ readily cancels from both sides of the equality. Similarly to 
section 4.1 we define <&ki = 4>ati,i + <&ki ,2 + ^ai ,3 where each is associated with 

the temperature gradient in direction i (as given by a right-handed set with x\ being 
the direction normal to the boundary) and is a solution to the Boltzmann-type equation 
(34) with boundary condition: 

dT G0 


— fV 


dxi 


+$/£!,i|n=0 = — 


-n 


j 9 T G o 


r/=0 


71 JQ'< 0 


dx; 


T )=0 


T ^Kl,i\i^=0 I l d fi , 


for fli > 0 
(46) 


We find that solutions (39) satisfy the above conditions for i = 2 and i = 3 respectively, 
and do not impose any condition over the tangential derivatives of T G q. For i = 1, (46) 
results in 


- ( 3+^! 


dT G0 


dxi 


= — $Kl,l|n=0 “ 2 


77=0 


[ ^i,i|r,=o^i^i, for > 0 (47) 

Jn\ <0 


The only solution possible with this boundary condition is <&a:i,i | ?;=0 = ( dT G o/dxi)\ v=0 = 
0. This can be seen by noting that if (dT G o/dxi) |^=o 7 ^ 0, multiplying the above equa¬ 
tion by Hi and integrating over 0 < fii < 1 yields ^ki,i|? 7 =o^i^i 7 = 0 , which 
is impossible (this can be seen by starting from the equation governing 4> 1.1 -of the 
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type (|34[) and integrating over 0 < r) < oo and —1 < fR < 1 and using the condition 
&ki,i(v —> oo) —> 0). We thus conclude that T G q must satisfy the boundary condition 


OTgo 

dn 


= 0 , 
Xf, 


(48) 


which is agrees with the Neumann boundary conditions associated with adiabatic bound¬ 
aries. 


5 Order 2 boundary layer analysis 

5.1 Order 2 analysis for prescribed temperature boundaries 


The second order correction &k 2 must be solution of (35) for i = 2, namely 
8&K2 


fb 


= (Kn) 


1 dr) 

with the boundary conditions 


£(®K2) — ®K2 


Kn 


- ^2 


d®Ki 

Bxo 


+ 


<94>xi\ 

dx 3 ) 


(49) 


$K2\rt=0 = — ®G2\r)=0 = ~ 1 g 2 | t )=0 + /t ^ x X] ^ G1 


Kn" 


77=0 ( Kn > 






d 2 T G0 

dxidxj 


r /=0 


(Kn) z —' ‘ dx 

' ' l . v, 

for fii > 0 (50) 

Here we note that the derivatives of the first order boundary layer which appear in the 
RHS of (49) introduce four volumetric source terms in the governing equation. 

The boundary condition (50) includes three terms with first order partial derivatives 
of Tq\ and nine terms with second order derivatives. Taking into account the four 
source terms on the RHS of (49), we introduce sixteen constants such that the order 2 

en 

d 2 T G0 


“temperature jump”, T G 2 \ v = 0 i may be written as 
3 3 


Tc'2 1 ?7=o = yx 


dT G1 


i=l 


dxi 


+ 22 g 8 
7 1 = ® i,j =1 


dx^Xj 


+ 22 ~ g 8 

r?=0 ij =2 


d z T G0 


dx^Xj 


(51) 


77=0 


We accordingly introduce sixteen boundary layer functions such that the total order 2 
boundary layer may be written as: 


< &K2 = 22 ^K2,i 


9T CI 


i =1 


dxi 


+ 22 

V=0 i,j =1 


d 2 T ( 


GO 


13 dxjdxj 


+ 22 ^ K2 

v=° i,j =2 


d 2 T ( 


GO 


13 dxidx-j 


(52) 


77=0 


The 16 unknown coefficients and boundary layer functions can be determined using a 
combination of numerical and analytical techniques; these are discussed in Appendix [B] 
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Here we summarize the final result, which, conveniently, is quite compact. The second 
order temperature jump is given by the condition 


-^G2177=0 = Cl 


dT G1 


dn 


(53) 


r /=0 


Due to its simplicity and compactness, this result lends itself particularly well to implicit 
application of boundary conditions; this is discussed in section |6.2| The analogy to the 
order one temperature jump extends to the temperature boundary layer that is given 
by 

rp. 9T G1 

J-K 2,1 = T K\, 1 


dn 


(54) 


77=0 


In addition to this temperature boundary layer, the analysis yields a second order heat 


flux boundary layer. It may be calculated analytically by inserting expression (52) for 
<&K2 into 

C ^K 2 V g d 2 ftdLO, (55) 


qjrafa) = [ 


An 


which can be written in terms of incomplete Gamma functions. Validation of these 
results can be found in 1221. 


5.2 Order 2 analysis of a diffusely reflective boundary 


In section 4.2, we resorted to an analysis of the order 1 boundary layers to obtain the 
order 0 boundary condition, and showed the latter amounts to the well-known Neumann 
boundary condition. Similarly, we here proceed with the order 2 analysis in order to 
find the boundary condition for the order 1 temperature field. 

Inserting (28) in (|44|) and introducing a boundary layer term yields, for Oi > 0 and 


for all frequency/polarization modes: 


T G 2 1X;, - Tjy ' ft ■ V x T G i| Xi) + 


(Kn> 

- f n\ (T G2 \ Xb - 

w Jn\< 0 V 


Kir 

W) 


; n ■ v x (n ■ v x t go ) | X6 + 4 >a'2|x 6 = 


K n K n ^ 

■ v x r G1 | X6 + -J^n' • v x (sr • v x t go ) 


(Kn 


(Kn) 




+ ^K2|x b ) d 2 n' 


Moving to the coordinate system (xi,X 2 ,xs) and the stretched coordinate 77 , we first 
note that in (56), the derivatives d 2 T G $/(dxidx\) | Q are zero for i = 2,3 because 
{dT GO /dxi)\ v= o = 0. 

Boundary layer <&K 2 may be decomposed into 4 components, 4> K 2 ,h ^A' 2,21 ^A' 2,3 
and $A 2 , 23 - Components &K 2;2 and <I>^ 2,3 are similar to the order 1 boundary layers 
*1^1,2 and (see expression (39)), with the only difference being that T G 0 is replaced 

by T G i. Component /V 2,23 corrects the anisotropic mismatch associated with the bulk 


(56) 
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term 2 Q J 2 ^ 38 2 T G q/{ dx 2 dx 3 ). It is a solution to the ID Boltzmann equation (34) with 
boundary condition 

<Sa- 2 , 2 3 |, =0 = -2ShSh (57) 

for Di > 0 , and 0 at infinity, and is therefore given by 

d 2 T G0 


r ]=0 


®K2,23 = — 2D 2 D 3 


8x 2 8x 3 


71=0 


** ' W 1 H(Ql) 


(58) 


Components &K 2 , 2 , ‘h K 2,3 and &K 2,23 do not contribute to a temperature jump or 
(temperature) corrective layer, but they do contribute to the heat flux boundary layer. 
The last component is solution to the following problem: 


n- 


9&K2, 

<9?? 




Kn 


i=2 


(Kn) 


dx 2 


exp 


r )=0 


-ry(Kn) \ 
fiiKn ) 




Kn 2 


T )=0 


+ \ 2 1^1“- 


(Kn) 


1\ 8 2 T< 


GO 


2) dx 'l 
2 o' 


+ 


7/=0 


Kn 

W) 


2 3 


ED 


i =2 


ff(fil 


1\ 3 2 T ( 


GO 


<9x\ 2 


+ |? 7 =o = —- / Hi < f > /<' 2 ,i|r)=o ^^ / 1 for Di > 0 and all w,p 

n Jn[<o 


lim $K2,i{to,u,P,v) = 0 

T )—^OO 


(59) 

Although we could solve problem (59) using the numerical method described in 122], 


we will here directly find the value of 7 without specifically calculating <&K 2 ,i- We first 


proceed by multiplying the boundary condition (second equation of problem (59)) by 
Di and integrating over the half sphere described by fK > 0 to obtain 


/ Hl < f > A'2.l|r;=od 2 fI = 

Jn 


47 t Kn OTqi 
3 (Kn) dx\ 


(60) 


rj =0 


We also multiply the first equation of problem (59) by V g C UiP and integrate it over all 
frequencies and solid angles and 0 < ij < 00 to obtain 

/ C U j t pVg^li < ^x2,i 17 )— 7oodojd 12 I C LJ: pVgQi<&K2,i\ri=odojd 12 

Jq,uj,p Jn,uj,p 


7 r 

'4 


f T T Kl12 ^ ,, 4-^00 ^ 

/^(Kn^^g 9x 2 (61) 


Since 4>a'2,i tends to 0 at infinity and V 2 Tgq = 0, we deduce the jump relation 


8T G1 


dx 1 


= 7 


9 2 T G q 


r /=0 


3s 2 


r )=0 


(62) 
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with 


3 L, P Kn2V 9 C co,pdu 
16 <Kn) f^ p KnV g C U)P dui ’ 

which can be rewritten in the form 

3 L, P V 9 T ‘ 2c ^P d “ 

' “ (A) faj,p VgTCu'pdlU ' 


(63) 


(64) 


In the single free path model, 7 = —3/16. Validation of this result can be found in 22 


Note also that the approach that we used for finding 7 may be used for finding the heat 
flux associated with the boundary layer $k2,i- 

A note on the physical interpretation of (62) At first glance, the boundary condi¬ 
tion (62) seems to suggest that energy is not conserved since the net heat flux into 
the (diffusely reflective) boundary is not zero. In fact, contrary to appearances, this 
form ensures energy conservation at the boundary. This can be seen by considering 
that d 2 Tco/dx\ 7 ^ 0 (only possible in two or three dimensions) implies variations in the 
temperature gradient along the boundary, which in turn implies variations in the heat 
flux along the boundary due to first-order kinetic boundary layers (see (39)). Impos¬ 
ing energy conservation at the boundary reveals that (62) exactly balances the terms 
resulting from gradients along the boundary 1221 . 


6 Summary and discussion of results 

We have derived the continuum equations and associated boundary conditions that 
provide solutions equivalent to those of the Boltzmann equation up to second-order in 
Knudsen number for steady problems. This derivation shows that the governing equa¬ 
tion in the bulk, up to at least second order in Knudsen number, is the steady heat 
conduction equation with the bulk thermal conductivity. Kinetic effects, always present 
at the boundaries due to the inhomogeneity introduced by the boundary and the con¬ 
comitant mismatch between the distribution introduced by the kinetic (Boltzmann) 
boundary condition and the distribution function in the bulk, become increasingly im¬ 
portant (can be observed in larger parts of the physical domain) as the Knudsen number 
increases. Fortunately, these kinetic effects can be systematically described and incorpo¬ 
rated into the continuum solution relatively straightforwardly via the addition of kinetic 
boundary layer functions that are universal for a given material and material-boundary 
interaction model. 

We have studied two types of kinetic boundary conditions: prescribed wall tem¬ 
perature and diffuse reflection. We now summarize the procedure for obtaining the 
temperature and heat flux fields for an arbitrary problem of interest. 

Prescribed wall temperature boundary condition : Let Tfc(xft) denote the prescribed 
temperature along the system boundary denoted by xj, with boundary normal n. Ac¬ 
cording to the asymptotic theory, the temperature and heat flux fields can be calculated 
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from 


T(x) = T 0 (x) + (Kn )(T G1 (x) + T K1 (x) ) + (Kn) 2 (T G2 (x) + T K2 {x) ) + 0((Kn) 3 ) 
q"(x) = (Kn)( q G1 (x) + q" a (x) ) + (Kn) 2 ( q' G2 (x) + c& 2 (x)) + 0((Kn) 3 ) 

where 


• T 0 (x) is solution to V 2 To = 0 subject to To | Xb = 7&| X6 

• T G i(x) is solution to V 2 T G i = 0 subject to T G i| X(j = | x& 

• T G 2 (x) is solution to V 2 T G2 = 0 subject to T G2 \ Xb = ci^^\ Xb 

• T k i(x) = T K i } 1 (ri)^\ Xb 

• T K 2 (x) = T, <{A {ri) d ]£ l \ Xh 


(Kn)q Gi = — kV x /T G i_i, i = 1,2 


<&i (*) = ZLf. 


u4 tt 


^a^ K1 4? ? )dujd 2 n §g 


ej with ^ki i,i = 2, 3 given 


by (39). 




q'A' 2 ( x ) = Iui, P ,n % E V g $A- 2 (? 7 )da;d 2 f 2 with $ K2 given by ((52]) 


We recall here that the coordinate r) is a stretched (by (Kn) -1 ) version of the local normal 
to the boundary. The boundary layer functions 7 x 1 , 1 ( 77 ), ^ Ki,i(v ) and T K2(d) are 
unique (universal) for each material and material-boundary interaction model. Figure 
[T] shows results for 7x1,1(77) for three material models. The method for calculating this 


function is described in detail in 22 . The boundary layer functions ^Ki,i(v) and ^a' 2 (t?) 


are known analytically. We also note that due to the absence of kinetic boundary layer 
corrections, at order zero T G o = Tq. 

Diffusely reflecting boundary: In the case of a diffusely reflecting boundary located 
at Xfr with normal vector n, the temperature and heat flux fields can be calculated from 


T(x) = T 0 (x) + (Kn)T G 1 (x) + O((Kn) 2 ) 

q"(x) = (Kn)(q G1 (x) + q^i(x)) + (Kn) 2 (q G 2 (x) + q^ 2 (x)) + 0((Kn) 3 ) 


where 

• To(x) is solution to V 2 To = 0 subject to ^j)| X( , = 0 

• T G i(x) is solution to V 2 T G i = 0 subject to with 7 given by 

©• 

• (Kn)q^ = -KV x /T Gi _i, * = 1,2 
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e ? ; with = 2,3 given by 


• q'xiW = EL L, P .n ( ^ir L n i yKi,idujd 2 rL 

©• 

• q^ 2 (x) = f^psi %f-Vg$K 2 (r/)dujd 2 n with the components of <f>A '2 given in sec¬ 
tion 15.21 

We note here that ^Ki,i is identical to the corresponding boundary layer function that 
appeared in the prescribed-temperature boundary condition case. We also note that 
due to the structure of the boundary-layer problem for the diffusely reflecting boundary, 
the first-order analysis yields a zeroth order boundary condition, while a second-order 
analysis yields a first order boundary condition; as a result the asymptotic solution for 
the temperature terminates at first order in (Kn). 

We see that, in both cases, the ’’traditional” Fourier description corresponds to the 
zeroth order solution. 


6.1 A one-dimensional example 

In this section we consider a simple ID problem as a means of illustrating the appli¬ 
cation of the asymptotic theory to problems of interest. We consider a silicon slab of 
thickness L confined between two boundaries at different prescribed temperatures. Us¬ 
ing dimensionless coordinates, the boundaries are located at x\ = — 1/2 and x\ = 1/2 
and have deviational temperatures Tr and Tr, respectively. 

We recall that under the asymptotic analysis, the temperature field is given by 

T(x i) = T 0 (xi) + (Kn)(T G 1 (xi) + T Kl { Xl )) + 0« Kn) 2 ) (65) 


The order 0 solution straightforwardly reads 

T 0 (xi) = Tl + Tr + (Tr - T l )x i ( 66 ) 

since it is the solution of the heat conduction equation subject to no-jump boundary 
conditions. Therefore, the boundary conditions for the order 1 field are 

T G1 (x 1 = ±1/2) = ± ci -^ = ±ci (T L - T r ) (67) 

OX 1 


which results in 


T G i(xi) = 2 ci (Tl — Tr)x i 


( 68 ) 


The boundary layer (Tr — Tl)tki,i((x\ ± l/2)/(Kn)) contributes to the solution near 
the boundary at x\ = — 1 / 2 , while the function (Tr — Tr)tr'i j i(( 1/2 — ^i)/(K n )) con¬ 
tributes close to the boundary at x\ = 1/2. The resulting solution correct to order 1 (eq 
(65)) is plotted in figure [3] for (Kn) = 0.1 in the single relaxation time model and com¬ 
pared to our benchmark (adjoint Monte Carlo |l 8 | ) result. The agreement is excellent; 
we note in particular that even though the boundary layer correction is small at this 
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Figure 3: Order 0 (dot-dashed line), order 1 (dashed line) and order 2 (plain line) 
solutions compared to the solution computed by highly resolved Monte Carlo simulation 
at (Kn) = 0.1. 


Knudsen number, the temperature jumps are considerable and are accurately captured 
by the asymptotic solution. In contrast, the zeroth order solution (which neglects the 
temperature jumps) is clearly inadequate. 

If desired, calculation of T(x i) to second order in (Kn) proceeds by solving the heat 
conduction equation for T G 2 subject to the second order boundary conditions. Applying 


(53) to this problem yields 


with the solution 


Tg 2 (x 1 = ±1/2) = =Fci G1 = ±2c\(T r - T l ) 

OX 1 


T G2 (*i) = 4 c\(T r - T l )x 1 


(69) 


(70) 


The order 2 solution including kinetic boundary layers is also shown in figure [3] and 
clearly exhibits improved accuracy with respect to the order 1 solution. In fact, in this 


particular problem where only first derivatives are non zero, the process by which (70) 


was derived can be repeated for all orders without knowledge of the higher order jump 
coefficients, leading to an asymptotic solution that is, in principle, correct to all orders. 
In other words, for n > 1: T Gn (x 1 ) = (—2) n c”(Tft — Ti)x\ 

Summing all orders (provided 2(Kn)ci < 1), we obtain: 


Tg(x 1 ) ~T L _ 1 


Tr-T l 


= - + 


xi 


1 + 2(Kn)ci 


(71) 


21 















The boundary layer corrections of all orders can also be obtained (and summed) using 
the same process. For example, for the boundary at xi = —1/2, we obtain 


Tr(x i) _ (Kn) / xi + 1/2 \ 

T r -T l ~ 1 + 2(Kn)c 1 T/a ’ 1 V (Kn) ) ' 


(72) 


The second boundary layer (at x\ = 1/2) is obtained in an analogous fashion. This 
solution is asymptotically accurate to all orders, meaning that the error converges to 
0 faster than any power of (Kn); for a discussion on the error associated with the 
asymptotic expansion see jTf> . 

Figure [4j compares the order 1, infinite order and “exact” (Monte Carlo) solution 
for (Kn) = 0.4. The infinite order solution is in very good agreement with the exact 
solution, while the order 1 solution is clearly inadequate at this Knudsen number. 


6.2 “Implicit” boundary conditions 


In the rarefied gas dynamics literature 17 jump boundary conditions are frequently 
imposed in an “implicit” fashion (in the sense that the unknown is on both sides of 
the equation, resulting to what is referred to in the mathematical literature as mixed 
boundary conditions) thus avoiding the “stagerred” solution procedure shown above 
where the governing equation needs to be solved for each order. For example, a set of 
boundary conditions up to second order given by 


Ko| X|) = T b 


rp , dT 0 




and 


T G2\ Xb = a 


dT G1 


dn 


+ /3 


d 2 Tn 


X b 


dn 2 


x b 


may be imposed by solving V 2 T(j = 0 subject to 

dT G 


Tg Ix 6 -T b = a(Kn) 


dn 


+ /3(Kn 


2 d 2 T G 


x b 


dn 2 


X b 


One can show that these two approaches are equivalent (to order (Kn) 2 ) 
TgU = (To + (Kn )Tgi + (Kn ) 2 T G2 + ...) 


I Xb 


(73) 

(74) 

(75) 

(76) 

by expanding 

(77) 


and similarly for dT G /dn | X() and substituting into (76). Equating terms of the same 
orders of (Kn) we obtain equations (73), (74) and (75), at order zero, one and two, 
respectively. 

Clearly the implicit form relies on the jump coefficients (a, /3 , etc) remaining the 
same at each order (e.g. in (74) and ([75])). If the above condition is satisfied, in 
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addition to requiring less solutions of the governing equation, the implicit form has one 
more advantage: provided higher order derivatives (not included in ([76])) do not appear 
at higher order, the solution will be correct to all orders, since it is easy to verify that 


(76) then implies that 


T Gn+ 2| X6 = a 


dTcn+l 


dn 


+ P 


d 2 T, 


Gn 


*6 


dn 2 


(78) 




for all n > 0. 

This property can be illustrated with the example of section |6T| where a = c\ and 
(3 = (1: solution ®> can be obtained directly by solving d 2 To/dx 2 = 0 subject to 


T g k -T b = Cl (Kn) 


8 T g 


dn 


(79) 


Xf, 


Although an infinite order solution is always welcome, we also need to keep in mind 
that some fortuity was involved in this problem in which all higher derivatives of the 
solution are zero. In the general case, given that (3 = 0, we expect the implicit condition 


(79) to provide solutions that are accurate at least to second order and at most up to 


order m — 1 where m denotes the order of derivative featuring a non-zero jump coef¬ 
ficient. We close by noting that the implicit approach sometimes results in boundary 
conditions which feature derivatives of the same order as the governing equation which 
may raise questions about the well-posedness of the mathematical problem. As a res¬ 
olution to this paradox, we recall that the derivation process followed here (sections [4] 
and [5]) produces the staggered forms of the general type (73)-(75), which do not present 
posedness problems. In other words, the implicit form is used merely for convenience 
and should be discarded if any mathematical/numerical issues arise. 


6.3 A two-dimensional example 

In this section we use a two-dimensional example to illustrate the application as well 
as convergence properties of the asymptotic solution theory. Specifically, we consider 
a slab of material that is infinite but subject to a periodic temperature variation in 
direction x \; the slab has thickness 2 L in the transverse direction, with the associated 
dimensionless coordinate X 2 defined such that X 2 = 0 describes the median plane of the 
slab. The material boundaries at X 2 = 1 and X 2 = —1 are at the prescribed (deviational) 
temperatures T w cos( 27 txi/ 3 ) and —T w cos(27rxi/3), respectively. The inset of Figure[5] 
shows a contour plot of the order 0 solution. 

In what follows, we construct the asymptotic solution of this problem up to 0((Kn) 2 ), 
both using the “order-by-order” approach and the implicit approach discussed in the 
previous section. We will then compare these solutions with MC simulation results, 
both visually along the line x\ = 0 but also very precisely at location {x\ = 0 ,X 2 = 1) 
to compare the order of convergence of the asymptotic solution with the theoretically 
expected one. 
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Figure 4: Order 1 solution (dashed line) and “infinite order” solution (solid line) 
compared to the solution computed by a finely resolved Monte Carlo simulation for 
(Kn) = 0.4. At this Knudsen number the boundary layer contribution is clearly visible 
(the solution is no longer a straight line). 



Figure 5: Zeroth order, Monte Carlo and implicit asymptotic solution for the tempera¬ 
ture along the line x\ = 0 in the two-dimensional example considered in section 6.3, for 
(Kn) = 0.1. The inset shows a contour plot of the order 0 solution. 
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The order 0 solution for the temperature held is given by 

sinh(pp) 


To(x i,x 2 ) = T w cos 


S inh (f) 


( 80 ) 


The order 1 bulk temperature held can be obtained by solving the Laplace equation 
with the boundary conditions: 


T g i(xi,x 2 = ±1) = Tci(Kn) 


dTo 

8 x 2 


{xi,x 2 = ± 1 ) 


(81) 


resulting in 


2 vr 


(Kn)T G i(xi,x 2 ) = —T w (Kn)ci— coth ( — 


n {2ir\ /27rxi\ sinh (Pp) 

coth ( T J cos (—j s J (¥) ' (82) 


The order 2 bulk temperature held is then obtained by solving the Laplace equation 
with the boundary conditions: 


T G2 (x 1 ,X 2 = ±1) = T c l (Kn) ~ 7 r— {x\,x 2 = ±1) 

OXo 


(83) 


leading to 


(Kn) 2 T G2 (xi,x 2 ) =T w (Kn) 2 c 2 


2 vr 


coth 


2 vr\ 

J 


cos 


/ 2ttx\ \ sinh (pp) 
\ 3 j sinh (p) 


(84) 


The solution is complete to second order once the boundary layer contributions are 
added. The order 1 and order 2 boundary layer correction terms in the vicinity of 
boundaries x 2 = ±1 are respectively given by: 


' 27r ( 2 i r\ 

(Kn)Tja(xi, x 2 ) = T?Vria,i((l =F a; 2 )/(Kn))(Kn) — coth — 

' 2 ( 85 ) 

(Kn) 2 T A - 2 (x'i,x 2 ) = ±T w ta'i,i(( 1 T x 2 )/(Kn))ci(Kn) 2 coth 


As explained in the previous section, a solution of a similar order can be achieved by 
directly looking for the solution of the Laplace equation T G with boundary conditions: 


dTr 

Tg(x i,x 2 = ±1) = ^ci(Kn)—^(xi,x 2 = ±1) 

0X2 


( 86 ) 


This is the case here because, as shown in section 5.1, second-order derivatives do not 
appear in the jump conditions or the temperature boundary layer. Applying these 
“implicit” boundary conditions, we obtain 


T g (xi,x 2 ) = 


1 + ci(Kn)p coth (p) 


cos 


/ 2:/rxp sinh (pp) 


3 J sinh ( 


2tt\ 
3 ) 


(87) 
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The kinetic boundary layer corrections in the vicinity of the boundaries at X 2 = ±1 
are given by = F T A'i,i((l T X 2 )/{Kn))(Kn)dTQ/dx 2 (x\ = 0,X2 = ±1). Evaluating the 
combined (bulk and boundary layer correction) solution at (xi = 0, X 2 = 1), we obtain 


T 


1 + ci(Kn)^ coth (^) 


( 88 ) 


This solution is compared to a highly-resolved MC simulation result in Fig. [5] for the 
case (Kn) = 0.1. The material model used is the single-relaxation-time model defined 
in section [4.1.1| The MC solution was obtained using the adjoint Monte Carlo method 
described in 1181 and will be denoted Tmc below. 

Figure [6] plots \T M c(xi = 0,x 2 = 1) - T asymp t ot ic(a;i = 0,x 2 = 1)| for 3 asymptotic 
solutions, namely, the first-order solution To + (Kn)(TGi +Tk 1 ), the second-order solu¬ 
tion Tq + (Kn)(T gi + Tki) + (Kn) 2 (TG 2 +TK 2 ), and the implicit solution (88). The figure 


shows that the implicit formulation leads to an order 2 solution overall which addition¬ 
ally features slightly improved accuracy compared to the “regular” order 2 solution. As 
explained in section 6.2, the solution would be “infinite” order if no higher order deriva¬ 


tive appeared in the jump boundary conditions. The third-order convergence observed 
for the implicit solution seems to suggest that a non-zero jump coefficient appears in 
front of the third-order derivative (m = 3). 



Figure 6: Convergence of asymptotic temperature solutions at (xi = 0, X 2 = 1) in the 


two-dimensional example considered in section 6.3 
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7 Extension to time-dependent problems 


Although the analysis presented here has so far been limited to steady problems, ex¬ 
tension to unsteady problems is relatively straightforward. In the held of rarefied gas 
dynamics the Hilbert expansion has been extended to time-dependent problems by 
Sone 101 and Takata [29,30 , who showed that, other than the additional time-derivative 
in the governing equation, time dependence does not introduce any new physics up to 
order 1 in (Kn). 

In this section we show that this is also true for phonon transport for the case of 
prescribed temperature boundaries by introducing the dimensionless time-dependent 
Boltzmann equation 

St^ + P-V x d>= ( 89) 

at Kn 

where t is a dimensionless time, defined by t = t'/to, where to is a characteristic time 
of variation and the Strouhal number is given by 


sw - st - 


(90) 


We analyze cases where (St) ~ (Kn), where the average Strouhal number, (St), follows 
from an analogous definition to that of (Kn) in ([6]). The condition (St) ~ (Kn) can be 
rewritten as to ~ L 2 /k ~ (r)/(Kn) 2 , which implies an assumption of diffusive scaling 
in time. 

Expanding the time dependent function as in Eq. (14) results in the same forms 
for orders 0 and 1 (equations (18) to ©)• Differences appear at order 2. Specifically, 
the form of the order 2 solution reads: 


Kn _ 

$G2 = £(®G2) — /v ' ^ • V x Tgi - 
(Kn) 


StKn dT ( 


+ 


Kn 2 


(Kn) 2 dt (Kn) 


-fi • V x (O • V x Tq) (91) 


Applying the solvability condition (25) results in 
StKn dT 0 


a 


Ld,P 


u, P ,n 47rr V(Kn) 2 dt (Kn 


Kn Kn 2 

+ • V x T G i - 7r7T7 ^ft • V x (n • V X T 0 ) ) d 2 fldu = 0 


(Kn) 


(92) 


which, after integration, yields the heat equation for the order 0 temperature field: 

3Tq K o m 

~w = ^ Vx,T °- 


C 


(93) 


Applying the solvability condition to the order 3 solution similarly yields the heat equa¬ 
tion for the order 1 temperature field. Although not strictly needed for our purpose 
here, we may solve for £(4> G 2) in Eq. (91) by writing: 


1 f CuJ ' p $ G2 d 2 Ctdu; 


Tr 1 9 — — 

C .L :P .n 4tt 


(94) 
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which, combined with (91), yields 


£ (‘f* G2 ) = Tq 2 + 


C'(Kn) 




C u , r StKndu,^ - 


a 


Kir 


CJ,p 


du\/lT 0 (95) 


' u),p 


which in the general case differs from Tq 2 - We note that £($^ 2 ) = 2g 2 holds in the 
case where the relaxation time does not depend on frequency and polarization. 


The order 0 boundary condition was obtained in section 4.1 by noticing that the 
order 0 distribution matches the distribution emitted by the boundary with no boundary 
layer correction. Introducing time dependence does not modify this result. Therefore 
the Dirichlet boundary condition Too = Tb remains unmodified at order 0 in the time- 
dependent case. At order 1, we showed that the jump boundary condition emerges 
from the analysis of the boundary layer correction required by the mismatch between 
the order 1 bulk distribution and the boundary emitted distribution. As before, time- 
dependence does not modify the form of the order 1 bulk distribution. Therefore, the 


order 1 jump condition (41) remains unmodified in the presence of time dependence. 


Similarly, the derivation of the order 0 condition for diffuse reflective walls results from 
an order 1 analysis. The Neumann condition (48) is unmodified. The order 2 boundary 


layer analysis presented in section 5.2 that yields condition (62) requires a modification 
since the relation V^To = 0 is replaced by the diffusion equation. In this work, we did 
not proceed to analyze in detail how the order 1 boundary condition for diffuse reflective 
walls is modified. 

This shows that the theory developed in this article may be applied to time-dependent 
problems (exhibiting diffusive scaling in time) up to order 1 in the presence of prescribed 
temperature boundaries, with the only change being that the Laplace equation is re¬ 


placed by the unsteady heat equation (93) 


7.1 Application to a transient problem 


To illustrate and briefly validate some of the conclusions of the previous section, we 
consider here a square particle heated to a uniform temperature of 301 K and placed 
in a thermal bath at 300 K, such that its boundary is well described by a prescribed 
temperature of T b = 300 K. We also assume that the Knudsen number is small such that 
we can calculate the temperature held inside the particle by solving the heat equation 


dT 

~dV 


= 


(96) 


with the first-order boundary conditions derived in this work. For convenience we use 

(97) 


the ’’implicit form” described in section 6.2 


dT 


T(x = x b ) - T b = ci(Kn) — 

In Figure 7l we show a measure of this temperature relaxation process, namely 


=0 


T(t) - T b 


for (Kn) = 0.1, where T(t) denotes the temperature at the center of the particle. The 


28 
















Figure 7: Temperature at the center of a square particle after initial heating. 


heat equation solution was obtained using a finite difference scheme. Here we note that 
the particle center is sufficiently far from the boundary that no kinetic boundary layer 
correction is required. The material model adopted here is that of silicon with a single 
relaxation time ( c\ ~ 1.13). 

This solution is compared with results obtained using the adjoint Monte Carlo 
method presented in Ref. jT8]. We also show the solution obtained from the (tradi¬ 
tional) heat equation with the Dirichlet boundary conditions T(x = Xb) = 300K. The 
figure shows that the asymptotic solution is in excellent agreement with the MC solution, 
while, as expected, the traditional approach (with Dirichlet boundary conditions)-which 
corresponds to the zeroth-order solution-significantly overpredicts the particle cooling 
rate. 


8 Application to interfaces between materials 

The theoretical and numerical considerations presented in this paper are quite general 
and can be extended to a variety of problems where boundaries introduce “size effects” 
by injecting inhonrogeneity into the problem. A classic example of such a problem is 
the interface between two materials: the presence of the interface results in a temper¬ 
ature jump, already shown in this work to be the signature of the kinetic correction 
required due to the inhomogeneity associated with the presence of a boundary. In this 
section we show how the asymptotic theory enables us to rigorously relate the Kapitza 
conductance to the kinetic properties of the interface (e.g. reflection/transmission coef¬ 
ficients) . Our aim here is not to conduct an exhaustive study but rather to demonstrate 
the applicability of the ideas presented earlier. As a result, we will focus on one specific 
transmission model and the single relaxation time model. We assume the following: 

- The interface separating the two media, denoted a and 6, is sharp (infinitely thin) 
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and planar. 

- When a phonon encounters the interface, it is either reflected or transmitted. In 
either case, its traveling direction is randomized while it keeps the same frequency 
and polarization. We denote the transmission probability from material a to 
material b by \ab, while p a b = 1 — Xab denotes the probability of reflection at the 
interface while traveling from a to b. Similarly, Xba and pb a denote the transmission 
and reflection probabilities for travel from b to a, respectively. 


In what follows, we will use r a , C U)P)a , V g , a and 75 , Vg.b to denote the relaxation 

time, frequency-dependent specific heat and magnitude of the (frequency and polariza¬ 
tion dependent) group velocity in materials a and b, respectively. As before and without 
loss of generality, let us align the interface with the X 2 — X 3 plane (at x\ = 0) and let the 
positive xi direction point from material a to material b. In this notation, the kinetic 
boundary condition associated with the interface is given by 


V„ 


g,b~ 


Coj, p,b 


a 


^~®t\ Xl = 0 + = Xab*i L=0-^ a 

4 JOi>0 47r 


Vg^iCp^l 


Vg.a 


a 


uj,p,a 


<L 


a lxi=0" 


IQl <0 


Xba$b\ xi = 


in !<o 

j2 


Pba^- b \ Xl =o + C ^v g ^ l d 2 n 


+ 


/Oi>0 


Pab® 


a lxi=0“ 




(98) 

where superscript “ + ” (resp. “ ”) refers to particles moving in the positive (resp. 
negative) x\ direction. 

The order 0 solution in each material phase is solution to the Laplace equation 
V^To — 0 with the condition 7bo,aUi=o- = Tco,b\xi=o+ = To|xi=o at the interface. 
Replacing < I ) a| : r 1 =o- and <J>j,| a . 1=0 + by Tq\ Xi= q in (98) and performing the integrations, 
we obtain: 


Lg,fcC , u,p,feT’o|xi=0 — XabTo\xi=oCuj,p,aVg,a “I" /5baTo|xi=oCoj,p,6l / g,fe 

Lg,aCaj,p,aTo|xi=0 = Xba^O\x\=oC u) ,p )b Vg,b T Pab^C)|xi=oCoj,p,aI / g,a 


(99) 


which implies 


0 — XabC{jo, P) aMg,a \baCui,p ) b^gfi- 


( 100 ) 


The principle of detailed balance guarantees that the above is true for all co,p. Note 
that the condition ?GO,aUi=o- = ^GO,b|xi=o+ = To|xi=o does not determine the value of 
cj > 0 | a . 1=0 = Tq|xi=o- The additional required condition is given by heat flux continuity: 


on 

dx\ 


xi=0~ 


= Kb 


dTo 

dx\ 


xi=0+ 


( 101 ) 
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Following the procedure of section 4.1 we find that the order 1 solutions 

Kn 

$l,a = T].a ~ \ & ' VTo a 

(Kn) 


$1,6 = Tl,b ~ 


Kn b 
(Kn) 


(102) 


• VT f 


0 ,b 


cannot satisfy condition (|98j) without the introduction of boundary layers. Here Kiij 
denotes V g ^Ti/L, while (Kn) is a “reference” Knudsen number calculated from the prop¬ 
erties of one of the two materials (results are independent of the chosen reference). 

We introduce two boundary layer functions 4 'Ka and 4 'Kb-, and two constants c a and 
Cb, anticipating temperature jumps from the order 0 at the interface of the form 


rp | _ dTo,a 

l,a\ Xl =Q~ ~ C a dna 


Tj 


dTr 


1 . & lii=0+ ° b 


0,6 


drib 


xi=0 _ 


11 = 0 + 


(103) 


Limiting our analysis to variations only in the x± direction, we insert the order 1 solution 
(boundary layer included) in condition (98), to obtain 


0,6 


xi=0+ 


V g ,b^f^ (—Kn+ c b (Kn) + %(Kn)) U 1=0+ ^ 

= J Xab Cul ^’ a Vg,a (-KiIqQ', - c fl (Kn) - fjf a (Kn)) U 1=0 - 

dT 0 .b 


n\ ( m\ 


xi=0 _ 


I Pba CuJ ^’ b V g ,b (-Kn btt[ + cj(Kn) + ^^ b (Kn)) \ Xl=0 + ^ 


xi=0+ 




y 9 ,a 


CuJ ' p,a (-Kn a fii - c a (Kn) - f^a(Kn)) U 1=0 - ^ T °’ a 


dx\ 


xi=0 


f Xba C yf L V g ,b (-Kn+ cj(Kn) + % 6 (Kn)) | X1=0+ ^ 


/ Pab ~ 


a 


uj,p,a 


V gA (-Kn a n[ - c a (Kn) - % a (Kn)) | X1=0 - 


dx\ 
dT 0 , a 


1 = 0 + 


n\dn\ 


dxi 


n\d.n\ 


a;i=0 _ 


We then solve this boundary layer problem numerically to obtain the condition 

dn a 


f 1 . 6 1\ ,a — CK a 


dx 


(Kn) 


(104) 


(105) 


xi=0 


with c = c a /n a + Cb/nb, describing the first-order temperature jump across the interface. 


The numerical procedure used is described in 22 
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8.1 Validation 


We test the asymptotic solution method outlined here on a simple one-dimensional 
problem with the following features: 


- The total length of the system is 2 L. The two materials are aluminum (—1 < 
x\ < 0, hence, material a) and silicon (0 < x\ < 1, hence, material b). Here, we 
emphasize that we perform this calculation to validate asymptotic theory describ¬ 
ing phonon transport across the interface. As a result, the aluminum model used 
here does not include electronic transport, which leads to kai = 27.7 W/mK. The 
choice of aluminum was motivated by the fact that this metal is frequently used 
as a transducer in transient thermoreflectance experiments 31-34 and thus a reli¬ 
able and well understood Monte Carlo simulation model - a priority for validation 
studies - exists (35 for this material. We use Zai and % to denote the range of 


frequencies of the two material dispersion relations, respectively. We also use a 
constant relaxation time model in each material; specifically, we take r a = 10 -11 s 
in A1 and = 4 x 10~ n s in Si. 


- A temperature difference of 1 K is applied across the system by imposing a pre¬ 
scribed temperature of 301 K at x\ = —1, while the boundary at x\ = 1 is 
maintained at 300 K. We note that the prescribed temperatures are used here to 
impose a temperature gradient onto the system. They are in no way linked to the 
interface model. 


- We define (Kn) as the ratio between the mean free path in the silicon phase and 
L. We choose L such that (Kn) = 0.1. 


The phonon transmissivities at the interface x\ = 0 are adapted from the model 
described in 26 , which given a “target” interface conductance G (as input), pre¬ 
dicts 


L 


Xab — 


<^ex A inx Si ,p 


Ccj.p.Al^o.Al 


+ 


i 


LgI A1 ,p Cj,p,A1 K,A1 fui£X Si ,p Cj.p.SiK,Si 


+ J_ 

^ 2 G 


(106) 


for frequencies in Zai H Zg; (0 otherwise). Coefficients Xba are deduced from the 
principle of detailed balance. 


Due to the one-dimensional nature of the problem studied here and the absence of 
higher than first-order derivatives of temperature in either material, an “infinite” order 
solution is possible: it can be obtained by solving the following system of four equations 
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in four unknowns (Tai(xi = — 1), Tai(xi = 0 ), T%\(x\ = 0 + ) and Tgj(xi = 1)) 
' 1 - T Ai (xi = -1) = CAi(Kn) (T A i(xi = 0") - Xai(xi = -1)) 

dT Si 


T S i{xi = 0+) - T A i(.ti = 0 ) = c(Kn)«si 


dx\ 


£ 1 = 0 + 


«A1 


dT Al 


dx\ 


dT f 


= «Si 


Si 


£ 1=0 


dx\ 


(107) 


£1=0+ 


. Tsi{xi = 1) = csi(Kn) (T Si (xi = 0 + ) - T Si (xi = 1)) 

We emphasize here that the temperature jump relations at x = ±1 (first and last lines 
in (107)) appear only because of the particular formulation used here for imposing the 
temperature gradient, namely using prescribed temperature boundaries far from the 
interface. Here, but also in general, the dynamics of the interface are solely described 
by the second and third lines of (107), namely heat flux continuity and the temperature 
jump across the interface. 

Our numerical results are shown in Figure [8} The figure compares the temperature 
profile obtained with the deviational Monte Carlo method 18,27 to the order 0, order 
1 and “infinite” order asymptotic solution. The order 1 solution provides significant 
improvement with respect to order 0. After adding the corresponding boundary layer 
functions we find that the infinite order solution agrees very well with the Monte Carlo 
result. Using this model, we obtain the actual conductance value G = 108 MWm _2 K _1 , 
which is very close to the “target” value 110 MWm _2 K _1 used as input to the model 
described in [26]. Perhaps more importantly, we note that the MC simulation also 
predicts a conductance value (obtained by extrapolating the bulk temperature profiles 
in order to calculate the temperature difference at the interface) of 108 MWm _2 K _1 , 
which is in perfect agreement with the (infinite order) asymptotic result. By comparison, 
the diffuse mismatch model predicts an interface conductance of G = 343 MWm _2 K _1 . 
This is consistent with the fact that the diffuse mismatch model results in an upper 
bound for the interface conductance 36 


We note that the “infinite” order solution may not be available in the general, 
higher-dimensional case. Related treatments of “connection” problems associated with 
different carriers have appeared in 37-40 . 


9 Final remarks 

We have presented an asymptotic solution of the Boltzmann equation in the small 
Knudsen-number limit. The resulting solution provides governing equations and bound¬ 
ary conditions that determine the continuum temperature and heat flux fields in arbi¬ 
trary three-dimensional geometries. Our results show that, for steady problems, the 
equation governing the bidk temperature field up to second order in the Knudsen num¬ 
ber is the steady heat conduction equation. We also show that, up to first order in 
the Knudsen number, the equation governing the bulk temperature field in transient 
problems is the transient heat conduction equation. 
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Figure 8: Temperature profile in a ID system with an Al/Si interface 


Although this result is expected (at least to first order in the Knudsen number) 
courtesy of traditional kinetic theory analysis |2,i8] (expanding the distribution function 
about the local equilibrium and giving no consideration to boundaries), the present 
work additionally derives the boundary conditions that complement this equation so 
that the resulting solutions of this system are rigorously consistent with solutions of 
the Boltzmann equation. In particular, the present work shows that the constitutive 
relation is only valid in the bulk, while a few mean free paths from the boundaries 
kinetic effects are always present. These effects not only modify the local constitutive 
relation (which is no longer of the Fourier-type), they also have a significant effect 
on the bulk solution by modifying the effective boundary condition subject to which 
the heat conduction equation is to be solved. These effective boundary conditions 
are derived for a variety of kinetic boundary conditions and shown to generally be of 
the jump type thus explaining the temperature jumps at the boundaries previously 
observed and remarked upon ||8l 9,41 . We note here that the jump conditions are 


universal (non-adjustable), while the jump coefficients and kinetic boundary layers are 
universal for a given material and material-interface interaction model; in other words, 
they are independent of system dimensionality and once calculated they can be used in 
any geometry of interest. Tabulated data for the various boundary layers derived not 
available in analytical form are available upon request. 

These results provide no evidence or justification for modifying the material con¬ 
stitutive relation (thermal conductivity) as a means of extending the applicability of 
the traditional continuum description to the transition regime-the underlying physics is 
considerably more complex. According to the asymptotic theory presented here, in the 
regime (Kn) < 1 (strictly speaking (Kn) <C 1) solutions consistent with the Boltzmann 
equation are obtained using a thermal conductivity that is equal to the bulk value; the 
modified (typically reduced) transport rate associated with size effects due to boundary 
presence is captured by the additional resistance introduced by the jump boundary con- 
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ditions as well as kinetic corrections that are to be linearly superposed to the final heat 
conduction result. On the other hand, by virtue of the expansion considered here, this 
work pertains to breakdown and extension of the classical Fourier description due to 
the inhomogeneity introduced by boundaries. As a result, it does not treat kinetic ef¬ 
fects appearing in a spatially homogeneous material such as those arising from temporal 
variations that are fast compared to the relaxation time, or spatial variations that have 
characteristic lengthscales that are on the order of, or smaller than, the phonon mean 


free path, that are also of interest to the scientific community 1331,42,43 . Ultimately, a 


theory that captures both classes of kinetic effects under a unified framework needs to 
be developed; we hope that this work is a step towards that goal. 

Our results are extensively validated using deviational Monte Carlo simulations of 
multidimensional problems. Studies in rarefied gas dynamics |17| show that second- 
order asymptotic formulations are reliable to engineering accuracy up to Kn ~ 0.4 and 
in some cases, depending on the problem simplicity, beyond. Our numerical validations 
support this finding. 

We note that the theory presented here assumes boundaries to be flat (no curvature). 
Curvature effects are expected to introduce additional terms in the effective boundary 
condition expressions and associated boundary layer corrections [l6 ! . This will be the 
subject of future work. 

Due to its ability to capture the inhomogeneity in the distribution function as¬ 
sociated with presence of boundaries, the present theory lends itself naturally to the 
description of the Kapitza resistance and temperature jump associated with the inter¬ 
face between two materials. We have shown that the asymptotic description produces 
results that are in excellent agreement with deviational Monte Carlo simulations. In 
other words, given transmission and reflection coefficients at the interface, the asymp¬ 
totic theory may be used to predict the Kapitza resistance without any assumption on 
the form of the distribution in the interface vicinity. 


A Derivation of the governing equation for the order 1 
and order 2 bulk temperature fields 

In this section, we show that Tq i and Tq 2 are solution to the Laplace equation. We 


start with the case of Tq\ . We apply the solvability condition (|25j) to <f>G2 to obtain: 

Kn 2 

C’o ,,pV g U • V x ( T G2 - -j^il ■ V x 4or + ■ 

J,P,P / 1 ' 


[ C U)P V q fl ■ V x ( Tq2 - • V x Tgi + Cl • V x {ft • V X T 0 )) dud 2 n = 0. 

■L.v.n \ \ Kn ) \ Kn ) J 


Integration over d 2 Ct removes terms containing odd powers of fyielding 


(108) 




Cu * Va 7iC) Y = 0 

' n ' i =1 


dx 2 


( 109 ) 
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from which we conclude that 


v'T G1 = 0 


( 110 ) 


To obtain the Laplace equation for Tq 2 , we apply (25) to <F G 3 . After carrying out 


the angular integration and cancelling terms containing odd powers of 12* we are left 
with 


i u ,p, n \ (Kn) * dx'l (Kn ) 2 dxidxjdx k dx, 


Thus, in order to show that the Laplace equation holds for T G 2 , we need to show that 

( 112 ) 



r r)^T 

/ V 0 d 2 n = 0 . 

n ijki OXiOXjdXkdxi 


Performing the angular integration, we obtain 


E = x E tAIL = Tvjvjro = 0 am 


47r 


d 4 T 0 _4vr 2 , 


'n 




' dxidxjdxkdxi 


5 dxfdx 2 - 5 

1,3 1 3 


as desired. 

It appears that this procedure can be applied to all higher order terms (T G 3 , T G 4 , 
etc.). 


B Determination of jump coefficients and boundary layer 


functions in eqs ( 51 ) and ( 52 ) 


Coefficients d{ and g t j (in (f5l]))and functions ^ K2,i and T K2.ij (in (52)) are determined 
by boundary value problems of the same form as the ones discussed in section 4.1 


satisfying equation (34). The problems that determine the coefficients gij include the 
source terms from that appear on the RHS of (49). In the interest of brevity, we 


only discuss the ones associated with the source term — d&Ki/dx 2 - The remaining two 
(associated with the term — d&Ki/dxz) may be deduced by analogy. 

For i = 2,3 coefficient c/ 2 i is solution to: 


p 


d^>R 2 


2 i 


1- 


(Kn) 


(C{^R2,2i) - ^K2,2i^j + ^y^Pj e x p 


dr] Kn V v 

’&K 2 , 2 i{tt,u,p, r i = 0) + g 2i = 0, for Pi > 0 
lim ^K 2 , 2 i(fi, u,p,r]) = 0 

77^-00 


(114) 


where H denotes the Heaviside function. Two results can be obtained immediately: 
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- Coefficients d\, d 2 and d% are solutions to the same problems as ci, C 2 and C 3 . 
Consequently, they are equal and their associated boundary layers are the same, 


provided To is replaced by Tqi in Eqs. (39), (41) and (42) 


- Coefficients gij and gij for i 7 ^ j are zero. For instance, it can be verified that 


^ K2,23 = ' 


Kn 77 f-r]( Kn) 

0 for Oi < 0 


for fli > 0 


(115) 


is a solution of (114). Solutions for all g,.j and gij,i 7 ^ j can be systematically 
obtained by solving the associated problem without the K 2 ,ij) term and then 
verifying that C('&K2,ij) = 0 

We are left with five undetermined coefficients, namely gu, g 22 , # 33 , <722 and < 733 . 
These can be determined using the numerical approach described in 1221 (suitably mod¬ 
ified in order to accommodate the volumetric source terms which appear in the mathe¬ 
matical formulation). Instead of following this approach, here we prove that 
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2—1 


fjii 


d 2 T 0 


dx ? 


d 2 T 0 


3 

+ dx 2 
V=0 i =2 UX i 


= 0 


(116) 


r]=0 


and that, therefore, the temperature jump associated with the second order derivative 
is zero, while the boundary layer, although not zero, integrates into a zero temperature. 


The remaining five coefficients in relation (116), can be determined by finding the 
function \h that satifies: 


3T (Kn) 


Hiv — 
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(t(§)-t)+^ 
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d 2 Tr 


®(n,uj,p,ri = 0 ) = - 
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cte? 


rj =0 

for Hi > 0 


(117) 


r?=0 


lim 4'(n,a;,p, 77 ) = 0 

77—^OO 


Let 4' = Ylk=i where for k = 1,..., 5 correspond, respectively, to the five bound¬ 
ary layer functions that are the counterparts of the five temperature jump terms in 
relation (116). We proceed with a strategy similar to the one used above, namely, 


solve for each 4^ individually, ignoring the contribution of and then evaluating 

T('Lfc). In the present case £(4^) / 0 but X)fc=i £('hfc) = 0; more details can be found 


m 


22 . This proves that yh, 4^ is the solution of (117) with the specified source terms 


and boundary conditions, and that the resulting boundary layer satisfies the boundary 


conditions without requiring a temperature jump correction, that is, relation (116) is 
proved. 
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